#ifndef _SOLUTION_H
#define _SOLUTION_H

#include "GlobalDefines.h"
#include "Excitation.h"
#include "Geometry.h"

class Solution
{
public:
	Solution(Excitation* ex):Ex(ex){};
	virtual ~Solution(){};
	virtual void ExtractMatrixS(EquivalentSurface* equal, Geometry* object, emx::zMatrix& S){};		
	void ExtractMatrixT(EquivalentSurface* geo_to, EquivalentSurface* geo_from, emx::zMatrix& T);
	void Fill_inc(EquivalentSurface* equal, emx::zVector& tmp);
	void Fill_sca(EquivalentSurface* equal, emx::zMatrix& S);
	void Inc_Initial(EquivalentSurface* equal, emx::zVector& tmp);
	void Transfer(EquivalentSurface* geo_to, EquivalentSurface* geo_from, emx::zMatrix& T);
	bool GetResults(EquivalentSurface** equal, int Num, emx::real_data thetas, emx::real_data thetae, int thetan, emx::real_data phis, emx::real_data phie, int phin, int polartype);
protected:
	emx::Complex3D L(emx::Real3D point, TriangleMesh* geo, int k);
	emx::Complex3D K(emx::Real3D point, TriangleMesh* geo, int k);
	const int GetNode3(TriangleMesh* geo, const int edgeIdx, const int flag);
	const emx::Real3D GetRho(TriangleMesh* geo, const int edgeIdx, const int flag);

	inline emx::dComplex G(emx::Real3D R)
	{
		emx::real_data RLength = sqrt(dot(R,R));
		return exp(-emx::J*Ex->Wave_Num*RLength)/RLength;
	}
	inline emx::Complex3D Nabla_G(emx::Real3D R)
	{
		emx::real_data RLength = sqrt(dot(R,R));
		emx::dComplex tmp = -exp(-emx::J*Ex->Wave_Num*RLength)
			*(1.0+emx::J*Ex->Wave_Num*RLength)/(RLength*RLength*RLength);
		return R*tmp;
	}
	emx::real_data Farfld(EquivalentSurface** equal, int Num, emx::real_data theta, emx::real_data phi, int PolarType);
protected:
	Excitation* Ex;
};

class MoM_Solution:public Solution
{
public:
	MoM_Solution(Excitation* ex):Solution(ex){};
	~MoM_Solution(){};
	virtual void ExtractMatrixS(EquivalentSurface* equal, Geometry* object, emx::zMatrix& S);		
private:
	void OutsideIn(Object* obj, EquivalentSurface* equal, emx::zMatrix& A);
	void FillZ(Object* obj, emx::dComplex* B, bool isColMajor);
	void InsideOut(EquivalentSurface* equal, Object* obj, emx::zMatrix& C);
};

//return the 3rd vertex's index
inline const int Solution::GetNode3(TriangleMesh* geo, const int edgeIdx, const int flag){
	if (flag == 0){
		int pat = geo->GetEdge(edgeIdx)->LeftFacet->Idx;
		return geo->GetPatch(pat)->Vertex(0)->Idx
			+geo->GetPatch(pat)->Vertex(1)->Idx
			+geo->GetPatch(pat)->Vertex(2)->Idx
			-geo->GetEdge(edgeIdx)->BeginVertex->Idx
			-geo->GetEdge(edgeIdx)->EndVertex->Idx;
	}
	else{
		int pat = geo->GetEdge(edgeIdx)->RightFacet->Idx;
		return geo->GetPatch(pat)->Vertex(0)->Idx
			+geo->GetPatch(pat)->Vertex(1)->Idx
			+geo->GetPatch(pat)->Vertex(2)->Idx
			-geo->GetEdge(edgeIdx)->BeginVertex->Idx
			-geo->GetEdge(edgeIdx)->EndVertex->Idx;
	}
}

//return the vector \rho
inline const emx::Real3D Solution::GetRho(TriangleMesh* geo, const int edgeIdx, const int flag)
{
	if (flag == 0)
	{
		int pat = geo->GetEdge(edgeIdx)->LeftFacet->Idx;
		return geo->GetPatch(pat)->Center-geo->GetNode(GetNode3(geo,edgeIdx,flag))->Position;
	}
	else
	{
		int pat = geo->GetEdge(edgeIdx)->RightFacet->Idx;
		return geo->GetNode(GetNode3(geo,edgeIdx,flag))->Position-geo->GetPatch(pat)->Center;
	}
}

#endif